function volumetric_locking_demo
clc
%close all
addpath('../');
addpath('../../');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% create mesh/read mesh
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
projDir = '.';
projName = '';
fn = [projDir '/' projName '/' ...
    'volumetric_locking_demo_linear.msh'];
mg = fem2.gid.GidMeshGenerator(fn);
mesh = mg.getMesh2d();
M = fem2.Model(mesh);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
meshUtil = fem2.MeshUtilities(mesh);
c1 = fem2.Constraint([0 1]); %fix x
c2 = fem2.Constraint([1 0]); %fix z
leftedge = meshUtil.seekNodesOnSurface(1,0,0);
botedge = meshUtil.seekNodesOnSurface(0,1,0);
mesh.setConstraint(c1,leftedge);
mesh.setConstraint(c2,botedge);
% plotMesh(mesh);
traction = [0.09801714	0.995184727;
            0.290284677	0.956940336;
            0.471396737	0.881921264;
            0.634393284	0.773010453;
            0.773010453	0.634393284;
            0.881921264	0.471396737;
            0.956940336	0.290284677;
            0.995184727	0.09801714];
tractionedges = [1 2;
            2 3;
            3 4;
            4 5;
            5 6;
            6 7;
            7 8;
            8 9];
for i=1:8
    ldx = math2.ConstantFunctionRnToR(traction(i,1));
    ldy = math2.ConstantFunctionRnToR(traction(i,2));
    n1 = mesh.getNode(tractionedges(i,1)-1);
    n2 = mesh.getNode(tractionedges(i,2)-1);
    e = mesh.addEdge(fem2.Edge([n1 n2]));
    el = fem2.EdgeLoad(e,[ldx ldy]);
    M.addLoad(el);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% element declaration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nf = mesh.countFaces();
mu = 1;
nu = 0.499;
E = 2*mu*(1+nu);
t = 1.0;
rho = 0.0;
% ss = fem2.StressState.PLANE_STRESS;
ss = fem2.enu.State.PLANE_STRAIN;
cm = fem2.StVenantKirchhoffMM(E, nu, t, rho, ss);
for i=1:nf
    e = fem2.StructuralElement(mesh.getFace(i-1), cm);
    %e.setIntegrationType(fem2.enu.IntegrationType.REDUCED_INTEGRATION);
    e.setIntegrationType(fem2.enu.IntegrationType.SELECTIVE_REDUCED_INTEGRATION);
    M.addElement(e);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solve
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mesh.setNumberOfFields(2);
M.preAnalysis();
n = mesh.countDofs();
kk = inf.jlinalg.Array2DMatrix(n,n);
M.assembleKt(kk);
r = zeros(1,n);
r = M.assembleR(r);
u = kk.getValues()\r;
M.setX(1,u);
for i=0:mesh.countNodes()-1
    mesh.getNode(i).getUHatLocal()
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% post process
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
plotDeformedMesh(mesh,1);
hold on
plotExactSolution(mesh,E,nu);

rmpath('../');
rmpath('../../');
end

function plotExactSolution(mesh,E,nu)
nn = mesh.countNodes();
nf = mesh.countFaces();
adj = getAdjacentMatrix(mesh);
P0 = mesh.getNode(0).getPosition();
P1 = mesh.getNode(nn-1).getPosition();
a = P0(2);
b = P1(1);
p = 1;
for i=1:nn
    X1 = mesh.getNode(i-1).getPosition();
    r1 = norm(X1);
    ur1 = (1+nu)*a^2*b^2*p/(E*(b^2-a^2))...
        *(1/r1+(1-2*nu)*r1/b^2);
    dX1 = [ur1*X1(1)/r1 ur1*X1(2)/r1];
    for j=1:nn
        X2 = mesh.getNode(j-1).getPosition();
        r2 = norm(X2);
        ur2 = (1+nu)*a^2*b^2*p/(E*(b^2-a^2))...
        *(1/r2+(1-2*nu)*r2/b^2);
        dX2 = [ur2*X2(1)/r2 ur2*X2(2)/r2];
        if (adj(i,j)==1 || adj(j,i)==1) && i~=j
            line([X1(1)+dX1(1) X2(1)+dX2(1)],...
                [X1(2)+dX1(2) X2(2)+dX2(2)],...
                'color','black');
        end
    end
end
end
